### CGQ 11 Reanalysis
### 4 December 2011

library(RItools)
source("code/seqblock1.R")
source("code/seqblock2k.R")
source("code/mahal.R")
source("code/balanceCheck.R")

source("code/cobgrequi11prep.R")

vars.to.evaluate <- c("MiniIsVtrWhite", "MiniIsVtrFem", "isVtrHmEngl", "Educ", "MiniVtrAge")
vn <- c("isWhite", "isFemale", "English at Home", "Education", "Age")

## Supplementary Figure 5:
out <- balance.check(dat = x.pr.sort, vars.to.eval = vars.to.evaluate, vars.type = c(rep("disc", 4), rep("cont", 1)), vars.names = vn, cex.pval.axis = .7, tr.var = "isTreated", plot.file.qq = "balCGQ.pdf", plot.file.pvals = "balpCGQ.pdf", xy.line.color = "grey")

## Sequentially block CGQ data (optional -- loaded below):
d2.p.CGQ <- array(NA, 100)
for(i in 1:100){
	print(i)
	## given dat X, do an SB.      
	sb1 <- seqblock1(id.vars = "myid", id.vals = x.pr.sort[1, "myid"], covar.vars = vars.to.evaluate, covar.vals = x.pr.sort[1, vars.to.evaluate], tr.names = c("0","1"), assg.prob.kfac = 5, file.name = "~/Desktop/cgqTest.RData", verbose = FALSE)
  
  for(n.idx in 2:nrow(x.pr.sort)){
  	sb2k <- seqblock2k(object = "~/Desktop/cgqTest.RData", id.vals = x.pr.sort[n.idx, "myid"],  covar.vals = x.pr.sort[n.idx, vars.to.evaluate], file.name = "~/Desktop/cgqTest.RData", verbose = FALSE)
  }

d2.p.CGQ[i] <- xBalance(I(as.numeric(Tr)) ~ MiniIsVtrWhite + MiniIsVtrFem + isVtrHmEngl + MiniVtrAge + Educ, data = sb2k$orig, report = c("chisquare.test"))$overall$p.value
}

load("data/d2pCGQ100.RData")

## Figure 7 (right, color):
plot((1:100)/101, as.numeric(sort(d2.p.CGQ)), xlim=c(0,1), ylim=c(0,1), col="blue", pch=18, ylab = bquote(d^2 ~ "p-values"), xlab="Sequential Reallocation", axes=FALSE, main = "Voter ID Training Experiment")
axis(1, at=c(1, 20, 40, 60, 80, 100)/101, labels=c(1, 20, 40, 60, 80, 100))
axis(2)
abline(h=out$d2.p, col="red")

## Figure 7 (right, black & white):
plot((1:100)/101, as.numeric(sort(d2.p.CGQ)), xlim=c(0,1), ylim=c(0,1), pch=18, ylab = bquote(d^2 ~ "p-values"), xlab="Sequential Reallocation", axes=FALSE, main = "Voter ID Training Experiment")
axis(1, at=c(1, 20, 40, 60, 80, 100)/101, labels=c(1, 20, 40, 60, 80, 100))
axis(2)
segments(x0 = 1/101, x1 = 100/101, y0 = out$d2.p, col = "grey", lwd = 2)
